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ABSTRACT 


Decompositions  of  the  plane  into  disjoint  components 
separated  by  airves  occur  frequently.  We  describe  a  package  of 
subroutines  which  provides  facilities  for  defining,  building  and 
modifying  such  decompositions  and  for  efficiently  solving  vari- 
ous point  and  area  location'"|:flrbblei!fs."'"'  " '""^ 

Beyond  the  point  that  the  specification  of  this  package  may 

be  useful  to  others,  we  reach  the  broader  conclusion  that  well 

::c:    03    ,s:;~s-:;i.-.^;:'   ...  ^;:. 

designed  data  structures  and  support  routines  allow  the  use  of 
more  conceptual  or  non-numencal  portions  <^f^  mathematics  in 
the  computational  process,  thereby  extending' greatly  the  poten- 
tial scope  of  the  use  of  computers  in  scientific  problem  solving. 
Ideas  from  conceptual  mathematics,  symbolic  computation  and 
computer  science  can  be  utilized  within  the  framework  of  scien- 
tific computing  and  have  an  important  role  to  play  m  that  area. 
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1.  INTRODUCTION 

Many  physical  systems  involve  several  materials  separated  by  interfaces. 
Recently  we  have  developed  a  library  of  subroutines  for  specifying  and  mani- 
pulating such  interfaces  in  the  plane.  In  this  paper  we  describe  the  capabili- 
ties of  this  package.  Applications  of  the  package  to  front  tracking  Ln  fluid 
flow  problems  have  been  described  elsewhere.^'2'^''*'^'^'^'8.9,l0,ll,l2 

Beyond  the  point  that  the  specification  of  this  package  may  be  useful  to 
others,  we  reach  the  broader  conclusion  that  well-designed  data  structures 
and  support  routines  allow  the  use  of  more  conceptual  or  non-numerical  por- 
tions of  mathematics  in  the  computational  process,  thereby  extending  greatly 
the  potential  scope  of  the  use  of  computers  in  scientific  problem  solving.  Put 
in  other  terms,  it  is  seen  that  ideas  from  conceptual  mathematics,  symbolic 
computation  and  computer  science  can  be  utilized  within  the  framework  of 

scientific  computing. 

.'Ci::."/;.   ic'?  br.£  ccr-ilr.-  ■. 
We  illustrate  this  point  by  the  u§c  .^hidi^has  bf*en  made  of  this  interface 

package  in  the  front  tracking  method  referenced  above.   There  are  three  basic 

difficulties  which  must  be  overcome' in  developing  a  front  tracking  method. 

They  are  (a)  intellectual  coherence,  to  control  the  proliferation  of  special 

cases,  (b)  adaptive  grids  fdi:' dynamically  evolving  lower  dimensional  fronts 

and  (c)  analytic  knowledge"  of  sollition  singularities.  The  package  described 

here  is  an  important  aspect  of  the  solution  of  problems  (a)  and  (b)  in  the 

front  tracking  program. 

Consider  a  set  of  points  (nodes)  in  the  plane,  joined  in  some  way  by 
non-intersecting  directed  curves.  We  will  use  the  word  interface  to  describe 
such  a  system.  In  another  terminology,  an  interface  is  a  planar  grap'i  each 
of  whose  edges  is  a  curve.  Given  an  interface,  the  plane  is  decomp'-'sed  into 
disjoint  components  whose  boundaries  are  composed  of  interface  curves. 
We  associate  an  integer  component  value  (or  color)  with  each  such  com- 
ponent, and  thus  each  curve  acquires  associated  left  and  right  component 
values.  We  do  not  require  that  every  component  have  a  distinct  component 
value.  Indeed  all  components  may  be  assigned  the  same  value  if  desired. 
However  the  requirement  that  curves  be  non-intersecting  (except  at  their 
end-points),  ensures  that  a  curve  separates  at  most  two  components.  An 
important  operation  is  one  that  determines  the  component  at  an  arbitrary 
point  j:,>'  in  the  plane. 


As  a  simple  example  a  circle  is  an  interface  containing  one  node  (any 
point  on  the  circle)  and  one  curve  (joining  the  node  to  itself  with  counter- 
clockwise orientation,  for  example).  A  circle  divides  the  plane  into  its  inte- 
rior and  exterior  regions,  which  we  might  label  with  integers  1  and  2  respec- 
tively. The  curve  then  has  1  and  2  as  its  associated  left  and  right  com- 
ponents.  For  an  example  of  a  more  complex  interface,  see  Figure  1. 

The  package  of  subroutines  we  have  developed  provides  a  mechanism  for 
defining  and  manipulating  planar  interfaces  of  arbitrary  complexity  in  an  effi- 
cient way.  Thus  addition  or  deletion  of  a  point  is  effectively  an  0(1)  opera- 
tion. Similarly,  the  operation  of  determining  the  component  number  con- 
taining a  point  x,y  is  effectively  0(1)  in  most  cases.  Considerable  internal 
pre-processing  of  data  is  involved  to  achieve  these  efficiencies.  The  routines 
provide  for  automatic  storage  allocation  and  error  handling.  In  addition,  the 
user  may  tailor  the  package  to  his  own  needs  by  adding  entries  to  the  basic 
data  structures  describing  interfaces  and  by  supplying  subroutines  which  are 
called  automatically  after  each  intdrfa^'operation. 

The  package  described  has  been  implemented  in  the  C  programming 
language  but  could  be  easily  written  in  a  similar  language,  including  FOR- 
TRAN. The  package  is  structured  as  a  subroutine  library  and  thus  can  be  used 
as  a  component  in  a  larger  application.  Two  applications  have  been:  a  shock 
tracking  code  for  hyperbolic  fluidtltiws,^'ahd  an  elliptic  partial  differential 
equation  solver  for  discontinuous  Coefficient  probleins.^^ 


2.   DATA  STRUCTURES 

The  basic  data  structures  to  be  defined  are  called  INTERFACE,  NODE, 
CURVE,  BOND,  POINT  and  COMPONENT.  Corresponding  routines  allo'v 
objects  of  these  types  to  be  created,  copied,  modified  and  deleted.  Consider- 
able redundancy  is  built  mto  these  definitions  in  order  to  make  the  package 
more  flexible  and  efficient.  In  appendix  A  we  present  a  complete  set  of  such 
data  structure  declarations  in  the  C  language.  Similar  declarations  would 
apply  in  Pascal,  PL/I  or  other  languages  with  user-defined  types.  In  practice 
the  library  routines  manipulate  only  pointers  to  these  data  structures.  In  the 
following  sub-sections,  we  discuss  these  basic  data  structures  in  more  detail. 


Appendix  B  presents  a  complete  list  of  the  operations  allowed  on  these  data 
structures  along  with  a  brief  description  of  their  function. 


2.1.  INTERFACE 

An  interface  is  regarded  as  a  topological  object  consisting  of  a  set  of 
nodes  joined  by  non-intersecting  oriented  curves.  These  concepts  are  embo- 
died in  the  definitions  of  the  mutually  recursive  data  structures  INTERFACE, 
NODE  and  CURVE.  The  INTERFACE  data  structure  is  a  structure  containing  a 
list  of  NODES  and  a  list  of  CURVES  which  join  those  NODES  pair-wise. 

2.2.  NODE 

A  NODE  is  a  structure  consistipg  of  4.POENT,  describing  its  location,  and 
two  lists  of  CURVES,  the  out^curyei^aa.d^  th^  in_curves,  corresponding  to  the 
curves  which  start  and  end  at  that  node. 

.  ,  :.  -.   _  -jar ^o~.o::?  l  zz  v  r 

2.3.  CURVE 

A  CURVE  is  a  structure  containiiig  start^^and  end  NODES  and  left  and  right 
COMPONENT  values.  This  correspond?. tp  ^iie  picture  described  earlier  of 
directed  non-intersecting  curves  which  join  pairs  of  nodes  and  sepairate  topo- 
logical components  of  the  plane.  To  represent  geometry  as  well  as  topology, 
we  consider  curves  which  are  composed  of  elementary  segments  and  we 
represent  curves  as  doubly-linked  lists  of  segment  data-structures,  which  are 
called  BONDS.  Thus  the  CURVE  structure  contains  first  and  last  BONT)S  as  well 
as  start  and  end  NODES. 


2.4.   BOND 

A  BOND  is  a  structure  consisting  of  two  POINTS,  its  start  and  its  end,  and 
pointers  to  two  other  BONDS,  next  and  previous.  Sets  of  BONDS  are  thus 
naturally  organized  into  doubly-linked  lists,  allowing  the  BONDS  associated 
with  a  CURVE  to  be  traversed  easily  In  either  direction.  For  BONT)S  that 
represent  linear  segments,  this  information  suffices,  but  for  higher-order  seg- 
ments (eg  parabolic  or  cubic)  further  parametric  data  is  stored  in  the  BOND 


data-structures. 


2.5.  POINT 

A  POINT  is  a  structure  containing  x  and  y  coordinates  describing  its  loca- 
tion in  the  plane. 


2.6.  COMPONENT 

A  COMPONENT  is  an  integer  specifying  a  connected  component  of  the 
plane.  Every  point  in  the  plane  which  is  not  on  some  CURVE  of  the  INTER- 
FACE has  a  unique  COMPONENT  associated  with  it. 

2.7.  User  Extensions  o^  Datft  Strj^^^s; 

To  facilitate  the  tailonn;>  of  ffiii^pScJcagerto  various  specific  problems, 
each  data  structure  includes  a  USER  data  structure,  which  may  be  vacuous. 
These  USER  structures  are  defiij©5Jk:hy:  ina^o  definitions  in  a  separate  file 
usertypes.h,  which  is  typically  diffcr^jiD^chjapplication.  For  example  in  a 
heat  conduction  problem,  the  entry  ^juc/^s-    '-ri;    \o    -jq. 

# define    USE^_P6lNT  float  temp; 

in  usertypes.h  forces  the  allocation  of  one  real  variable  associated  with  each 
POINT,  to  represent  the  temperature  at  that  point.  User-supplied  routines 
may  be  supplied  to  manipulate  this  extra  variable,  see  below. 


3.   OPERATIONS  AND  USAGE 

Basic  operations  are  provided  to  create,  modify,  delete,  print  and  read 
objects  of  each  of  the  above  types.  To  illustrate  the  nature  and  use  of  these 
routines,  we  describe  in  detail  the  routines  which  support  the  highest  level 
data  structure,  INTERFACE.  The  support  routines  for  the  other  data  types  are 
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similax  and  are  listed  in  appendix  B. 


3.1.   Interface  Operations 

The  basic  n^rrERFACE  operation  is  makejnterface  which  allocates  storage 
for  an  INTERFACE  and  returns  a  pointer  to  its  location.  All  pointers  to  struc- 
tures in  this  INTERFACE  are  initialized  to  default  values  of  NULL.  A  record 
of  the  INTERFACE  is  stored  in  an  internal  linked  list  of  currently  active 
INTERFACES  and  finally  the  INTERFACE  is  recorded  as  the  current  INTER- 
FACE -  all  operations  that  modify  an  interface  act  only  on  the  current  INTER- 
FACE. 

The  routine  copyjnterface  makes  a  copy  of  an  interface,  including  a 
copy  of  all  its  CURVES,  NODES,  BONDS  and  POINTS,  and  sets  the  current 
INTERFACE  to  the  new  copy.  The  standard  way  to  delete  an  INTERFACE  is  to 
call  delete Jnterface,  which  frees  all  "^or^e"  used  by  an  INTERFACE  and  its 
substructures  and  removes  the  JNTERFACE  from  the  list  of  currently  known 
INTERFACES.  ^   ,^xi^J^3^  litst;  KIIEJ  i   -.  : 

Interactive  programs  which'ri€e{3  to  inf^  Interfaces  may  use  the  routine 
re  ad  Jnterface.  This  creates  an  INTERF-AGE  by  a  call  to  makejnterface  and 
then  prompts  for  input  of  the  various  CUR^'ES  and  NODES.  Similarly 
print  Jnterface  prints  an  INTERFACE  (and  its  NODES  and  CURVES)  in  a 
prescribed  format  -  both  binary  and  ascii  representations  are  available.  Con- 
versely, read_print Jnterface  is  the  inverse  to  this  routine  and  will  read  from 
a  file  a  formatted  INTERFACE,  printed  previously  using  print  Jnterface.  It 
calls  makejnterface  to  allocate  storage  for  the  new  INTERFACE,  places  the 
data  as  read  from  the  file  into  this  new  INTERFACE,  and  returns  a  pointer  to 
the  INTERFACE. 

To  control  which  INTERFACE  is  the  current  INTERFACE,  three  routines 
are  provided:  currentjnterface  returns  a  pointer  to  the  current  INTERF.ACE, 
set_current Jnterface  switches  the  current  INTERFACE  to  another  specified 
INTERFACE  and  exists  Jnterface  checks  that  a  given  pointer  is  really  a  pointer 
to  a  currently  active  INTERFACE. 


3.2.  Operations  for  Other  Data  Types 

There  are  similar  support  routines  for  the  data  structures  NODE,  CURVE, 
BO^a^  and  POINT.  To  allow  convenient  looping  over  POINTS,  BONDS,  CURVES 
and  NODES,  the  routines  next_point,  next_J>ond,  next_curve  and  next_node  are 
provided.  These  remember  che  last  such  element  that  was  processed,  and 
resume  from  there  at  the  next  call  or  may  be  re-initialized  to  run  through  the 
interface  again. 

For  a  complete  list  of  the  operations  on  the  various  data-structures,  see 
appendix  B,  which  also  gives  a  brief  description  of  their  functions.  All  of 
these  routines  do  considerable  error  checking.  In  addition  each  routine 
accomplishes  as  much  related  work  as  possible.  Thus  for  example 
makejcurve  not  only  allocates  a  CURVE  structure,  but  also  inserts  the 
CURVE  into  the  current  INTERFACE  CURVE  list,  into  the  injcurves  and 
outjcurves  of  each  of  its  end  NODES  and  initializes  the  first  and  last  BONDS  on 
the  CURVE.  Generally,  each  routinecah 'be' assumed  to  do  all  of  the  obvious 
book-keeping  automatically.  . -....ll-       , 

:::;    2:r\nr:'    ,2noil£3:!ibG;:r. 

3.3.  Surgery  on  CnrTCS 

Routines  to  modify  a  previously  created  interface  are  very  important. 
insert _pointJn_hond  will  split  a  BOND  into  two  BONDS,  by  insertion  of  an 
intermediate  POINT.  The  inverse  to  this  function  is  delete _start_of_bond, 
which  removes  a  POINT  and  replaces  two  adjacent  BONDS  by  a  single  BOND. 
Because  BONDS  are  stored  as  a  doubly  linked  list,  these  curve  modifications 
can  be  performed  in  0(1)  time.  At  a  higher  level,  there  is  another  pair  of 
inverse  functions,  splitjcurve  and  join_curve.  The  first  splits  a  CURVE  into 
two  CURVES  by  the  introduction  of  a  common  NODE.  The  second  joins  two 
CURVES  to  form  a  single  CURVE,  by  the  coalescence  of  the  NODES  at  the  end 
of  the  first  CURVE  and  the  beginning  of  the  second  to  form  a  single  NODE. 


3.4.   Superimposed  Boundaries 

Frequently  in  applications  a  rectangle  surrounds  a  region  of  interest. 
Given  such  a  rectangle,  the  routine  setjjoundary  finds  and  adds  a  set  of 
boundary  CURVES  to  an  INTERFACE  in  such  a  way  that  the  resulting  INTER- 
FACE decomposes  the  plane  into  connected  components,  with  one  component 
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being  the  exterior  of  the  rectangle.  The  boundary  of  the  rectangle  is  then  the 
union  of  the  added  boundary  CURVES.  The  routine  also  checks  the  validity 
of  the  COMPONENT  values  for  each  CURVE  which  meets  the  boundary  rectan- 
gle, i.e.  it  determines  whether  the  COMPONENT  assignments  are  mutually 
consistent. 

The  set_boundary  routine  assumes  that  all  POINTS  of  the  INTERFACE  are 
on  or  inside  the  rectangle.  A  related  routine,  zoom  Jnterf ace,  will  zoom  and 
clip  a  given  INTERFACE,  deleting  all  POINTS  which  lie  outside  of  some  given, 
possibly  rotated,  boundary  rectangle.  This  is  of  especial  use  in  various 
refinement  situations  where  a  sub-rectangle  of  a  computational  domain  is  sin- 
gled out  for  more  detailed  examination. 


3.5.   Carve  Intersectiona 


In  order  to  preserve  the  iijtegrity.  .pf  ^COMPONENT  assignments,  it  is 
essential  that  curves  of  an  interface  intersect  only  at  their  end-points.  As  a 
result  of  dynamical  modifications,  curves  may  intersect,  or  become  self- 
intersecting.  To  check  for  these  occurrences  a  routine  intersections  is  pro- 
vided which  searches  for  all  intersections  of  the  bonds  of  an  interface.  Each 
such  intersection  is  recorded  in  a  special  CROSS  data  structure  and  a  list  of 
these  crossings  is  returned  'hy'iniersectwns.  A  crude  implementation  of  this 
routine  would  require  6Xn-)  operations  Where  n  is  the  total  number  of  BONDS 
in  the  INTERFACE.  The  use  of  precomputed  and  stored  data  reduces  the  time 
for  this  routine  to  a  value  which  is  typically  0{n).  The  worst  case  speed  for 
this  routine  is  0{n^),  a  value  which  rarely  occurs  in  practice. 


3.6.    User  Extensions  of  Operations 

In  order  to  make  the  package  more  fle.Tibie,  each  routine  finishes  by  exe- 
cuting a  user-supplied  code  segment,  defined  by  a  macro  Jefinition  in  a  file 
called  userint.h.  Typically  this  code  will  operate  on  the  user-supplied  data- 
structures  mentioned  previously.  This  code  may  of  course  be  vacuous. 
Consider  the  previous  example,  where  each  POINT  structure  contains  an  extra 
USER^POENT  variable  defined  in  usertpes.h,  representing  the  temperature  at 
that  point.  Suppose  one  wanted  to  initialize  the  temperature  distribution  as 
new  POINTS  are  created  to  be  the  distance  of  the  point  from  the  origin.    This 
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would  be  accomplished  by  providing  in  file  userint.h  the  following 
USER_Point  addition  to  the  standard  PointQ  routine  used  to  create  new 
POINTS: 

^def.ne    USERJ'oint         p^temp  =  sqrt{sqr{p-x)  + sqr{p^y))\   . 


4.  THE  STORAGE  ALLOCATION  SCHEME 

The  routines  that  create  and  modL<'y  INTERFACES  all  use  a  common 
storage  allocator.  There  are  no  re£^.■ictions  an>'where  in  the  code  on  the 
number  or  sizes  of  INTERFACES,  NODES,  CURVES,  BONDS  or  POINTS.  If  the 
storage  allocator  runs  out  of  space  at  soni?  point,  the  routine  that  called  it 
will  return  an  error  value,  generally: 0  or. a..''rEj7uI  pofater  as  appropriate.  The 
low  level  storage  allocator  is  not  eaHf^^iicectly*  .The  I>rrERFACE  code  main- 
tains a  higher  level  storage  rlloij^^ox  wj^jck^dtepenses  space  as  required  and 
only  occasionally  calls  the  low  level  Qnesaoqnu  bri.: 

The  high  level  storage  scJieme  fesoeiaPtei  a  chain  of  large  blocks  of  space 
with  each  interface.  When  the  blocks'^or  the  current  interface  are  exhausted 
and  more  space  is  needed,  the  low  level  allccatcr  is  called  to  allocate  a  new 
large  block  (Chunk)  of  size  CHUNK.  The  Chunks  for  a  given  interface  are 
stored  as  a  linked  list.  This  procedure  tends  to  keep  each  interface  stored 
contiguously  and  prevents  unnecessary  fragmentation  in  the  low  level  storage 
allocator. 

For  internal  purposes  the  code  maintains  a  list  of  INTERFACE  TABLES. 
Each  TABLE  summarizes  the  storage  statistics  of  an  INTERFACE,  and  in  fact 
contains  the  INTERFACE  structure.  This  list  of  TABLES  Is  itself  maintained  as 
a  linked  list,  and  as  new  INTERFACES  are  created  space  for  their  TABLES  is 
first  allocated  from  the  low-level  storage  allocator.  The  code  is  structured  in 
such  a  way  that  none  of  the  user-visible  routines  ever  explicitly  reference 
these  TABLES. 


5.   COMPONENTS  AND  TOPOLOGY 

An  interface  divides  the  plane  Into  connected  components.  The  data 
structure  COMPONENT  defines  a  coloring  of  these  components  into 
equivalence  classes.  There  are  two  main  routines  which  locate  a  point  with 
coordinates  x,y  relative  to  the  interface.  The  routine  component  determines 
the  COMPONENT  in  which  x,y  is  located.  The  routine  nearest_interface_point 
finds  the  closest  BOND  and  CURVE  to  x,y  and  the  location  on  the  BOND  which 
is  closest  lo  x,y. 

Both  component  and  nearest_interface_point  are  called  frequently,  and 
their  efficiency  is  of  considerable  importance.  The  efficiency  is  achieved  by  a 
precomputation.  These  routines  then  .eiccess  stored  data,  and  typically  require 
0(1)  time.  .      .  .'.OV:   ■ 

The  idea  behind  the  pre-pTO<»sSing''is"to  introduce  a  superimposed  rec- 
tangular grid.  Normally  this  is  ddne^y-tii;e''uselp,  through  a  call  to  the  routine 
setjnterface_grid.  If  not,  the  cG.dte JChbdses^-ai  grid  itself.  It  selects  the  smal- 
lest rectangle  oriented  parallel: -t*''th(B  fi^orSiQite  axes  which  contains  all 
POINTS  of  the  INTERFACE  and  imposes'*  gri'd  on  that  rectangle.  The  pre- 
processed  stored  data  cqnsjsjS;  of  lis!lft,iafi]^NDS  and  CURVES  which  enter  each 
mesh  block  of  this  grid.  If-tbero  iar€^  Pq>BONDS  for  a  given  block,  then  the 
COMPONENT  is  uniquely  defined,  ^d;  isvalso  stored  with  the  precomputed 
data.  This  same  data  is  used  by  the  routine  intersections  for  efficiency  rea- 
sons. 

While  this  particular  pre-processing  mechanism  appears  adequate  for 
many  problems,  it  is  not  optimal,  either  in  computation  time  or  storage 
space.  It  is  not  difficult  to  replace  the  current  code  for  component  with  other 
pre-processing  schemes.  A  crude  point  location  algorithm  is  easily  devised 
based  on  comparison  of  location  relative  to  all  n  BONDS  of  an  INTERFACE, 
and  requires  0{n)  computation  time,  but  no  pre-processing  or  extra  storage. 
The  algorithm  described  here  consumes  0{n)  +  0{g~)  storage  and  pre- 
processing time  when  the  superimposed  grid  has  dimensions  g  by  g.  How- 
ever the  point  location  computation  cost  is  then  effectively  0(1)  in  n  for  typi- 
cal interfaces  we  have  encountered,  although  the  worst  case  computation  time 
is  still  0{n)  -  corresponding  to  the  case  where  the  complete  interface  is  con- 
tained in  a  single  grid  rectangle.  For  a  review  of  the  point  location  problem 
and  related  areas  in  computational  geometry  we  refer  to  several  recent  papers 
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in  the  computer  science  literature.^-''^'*'^^'^^  Future  versions  of  this  software 
may  well  incorporate  other  approaches  to  point  location.  It  is  important  to 
stress  that  what  is  required  is  algorithms  with  good  behavior  on  normal  inter- 
faces, even  if  their  worst-case  behavior  is  poor  (eg  close  to  linear). 


6.   INTERFACES  IN  HIGHER  SPACE  DIMENSIONS 

Based  on  our  experience  with  planar  interfaces,  we  make  here  several 
comments  on  the  problem  of  defining  data  structures  and  support  routines 
for  interfaces  in  higher  space  dimensions.  We  adopt  the  point  of  view  of 
computational  continuum  mechanics,  where  the  interfaces  designate  material 
boundaries,  phase  boundaries  or  discontinuous  waves,  but  there  are  a  number 
of  other  and  unrelated  applications  such  as  graphics  for  these  ideas.  An  inter- 
face in  R"  will  consist  of  a  number  of  smooth  surfaces  and  their  boundaries, 
boundaries  of  boundaries  etc.  The  smooth  surfaces  will  in  general  have  an 
arbitrary  topology,  so  that  they  cannot  be  parameterized  by  an  open  subset  of 
R^  or  in  other  terminology,  by  a  sfe^fii^  cobrdiiiate  chart.  The  surfaces  of 
dimension  j  will  in  general  share  s6me  ""common  boundaries  of  dimension 
;'  -  1.  Thus  the  data  structure  for  a  numerical  description  of  an  interface 
will  consist  of  a  list  of  smooth  surfaces  of  various  diincnsions  together  with  a 
list  of  boundary  relations  among  these  surfaces.  Each  surface  of  dimension  j 
should  include  a  pointer  to  the  surfaces  of  dimension  j  +  I  which  it  bounds 
and  to  the  surfaces  of  dimension  y  —  1  in  its  boundary  together  with  the  rela- 
tive orientation  of  these  surfaces.  We  note  that  the  boundary  of  an  interface 
of  highest  dimension  j  is  an  interface  of  highest  dimension  j  -  1,  so  that  for 
the  case  n  =  3,  ;  =  2,  the  boundaries  are  given  by  the  data  structures  defined 
in  this  paper,  with  the  simple  addition  of  an  extra  real  variable  z  to  the  POEnT 
data  structure. 

There  are  three  requirements  which  determine  the  data  structure  for  the 
description  of  a  single  smooth  surface  in  an  interface.  A  small  perturbation  of 
the  points  (such  as  would  occur  in  a  single  time  step  of  a  time  dependent  evo- 
lution) should  not  require  a  rrparameterization  of  the  surface.  Reparameteri- 
zation,  when  required,  should  be  a  local  operation  and  the  local  reparameteri- 
zation  should  be  accomplished  in  0(1)  time.  The  grid  on  the  surface  should 
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bc  suitable  for  use  in  a  finite  difference  or  finite  element  computation  and  in 
particular  should  not  have  a  bad  aspect  ratio.  Because  the  surface  may  have 
an  arbitrary  topology,  there  will  in  general  not  be  a  global  rectangular  pattern 
to  the  index  structure  of  the  surface  elements.  Also  because  of  the  require- 
ment of  invariance  under  small  perturbations,  the  surface  elements  will  not 
be  regular  polygons.  Thus  it  is  natural  to  suppose  that  the  the  surface  ele- 
ments are  triangles  (simplices,  for  j  >  2,  n  >  3).  The  aspect  ratio  require- 
ment means  that  the  neighboring  triangles  do  not  differ  too  greatly  in  size 
and  that  each  individual  triangle  does  not  have  an  extreme  aspect  ratio.  It 
follows  that  the  number  of  triangles  meeting  at  each  vertex  should  be 
bounded  by  some  small  number.  Finally  in  order  to  support  stencil  operators 
of  a  finite  difference  scheme,  it  is  desirable  for  each  triangle  to  know  which 
triangles  are  its  contiguous  neighbors.' 

r.  :--;^  :.i  L\:m  v. 
...■      '■   lajoorcv;  \'>  -itr.'.:- 
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APPENDIX  A:   INTERFACE  DATA  STRUCTURES 

We  present  here  a  complete  set  of  interface  data-structures  in  the  C 
language.     These  structures  implement  the  data  types  described  in  section  2. 
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In  the  simplest  case,  each  of  the  USER  entries  in  these  structures  can  be 
defined  to  be  vacuous  by  adding  definitions  such  as: 

^define    USERJ>OINT    ; 

These  structures  are  mutually  self -referential.  They  are  first  defined  as 
structures  with  names  beginning  with  an  underscore,  and  the  actual  type 
definition  occurs  subsequently  in  a   C  ryp^cftf/ statement. 

typedef  int  cONffONENT; 


struct  .POINT  {  ■  ■  ■   ■'   ~.-Vl-.' 
float  x; 

Hoat  y;  ,..,    dbGL!?M    It. 

USER.POINT  .i:'J'.  .rj  .lov  ,t-:- 


}; 

typedef  struct  .POLVT  POLVT; 


._i  ifMirriio  -ganoic 


.'-.■'    "Ad':  :.  l>s- 
:;  '1 0    '  -  •'.  2  0  k-i'.'-'.  ;■,•  c  xi  c 


strua  30ND  { 

struct  _BONXi   'next; 
struct  JONT)    *prev; 
POLNT   'start; 
POINT  'end; 
float  length; 

USERJOND 

}; 

typedef  struct  _BOND  BOND; 
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struct  .CURVE  { 

strurt  INTERFACE  •interface; 
struct  J^ODE    'start; 
struct  _NODE    'end; 
COMPONENT    left_comp; 
COMPONENT    right_comp; 
int  boundary; 
BONT)       'first; 
BONT)     'last; 
int  num^oints; 
USER_CUR\'E 

}; 


struct  ^\ODE  {  ,, ,  ^^^  ,,_^  ^.^^  ,,  ^_ 

Struct  INTERFACE  •interfaceL„,.  ^..-.  l,_.  »,.^',^- 
POINT      'posn;  -d^.^c—^Z  r^  iJirqiJ 

struct  .CURVE    **in_airves;     -jj^;— ^^rj  ^  i^yjrr.uC, 
struct  _CUR\'E   ••out.curves;   -,^  x  jir^L'c  -.b:'" 


int  boundary; 

USERJ^'ODE 


Struct  INTERFACE  { 

struct  js'ODE     "nodes; 
Struct  _CL'R%'E   "curves; 
/•  Internal  Variables:  */ 
struct  Table  'table; 
int  modified; 
USER.INTERFACE 

}; 


typedef  struct  _CURVT  CL'RVE; 
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typedef  jtnict  _NODE  NODE; 

typedef  struct  .interface  INTERFACE; 


APPENDIX  B:   INTERFACE  OPERATIONS 

We  provide  here  a  complete  list  of  available  operations  on  C^TTERFACES 
and  their  sub-structures.  Most  of  these  functions  rerum  an  error  status  to 
indicate  incorrect  arguments  or  usage. 

Operations  on  Interfaces: 


makeJnterfaceO 
copyjnterfacei) 
delete  Jnterfacei) 
read_interface{) 
print  JnterfaceO 
read_print  Jnterfacei) 
currcntjnterfacei ) 
set_current_interface( ) 


Allocate  an  INTERFACE  structure 
Make  a  copy  of  an  INTERFACE 
Delete  and  free  stora'ge  fc3r  INTERFACE 
Inputs  an  INTERFACE  interactively 
Outputs  a  Formatted  INTERFACE 
Reads  output  of  print_interfaceQ 
Reports  the  current  INTERFaCZ 
Changes  the  current  INTERFACE 


make_node() 

copy_node() 

delete_node() 

print  _node() 

is_boundary_node( ) 

Ls_node_of_closed_cur\'e_only( ) 


Allocates  a  NODE,  adds  to  INTERFACE 
Copies  a  NODE  to  current  INTERFACE 
Deletes  a  NODE  from  INTERFACE 
Outputs  a  Formaned  NODE 
Checks  if  Node  on  Boundary 
Checks  if  a  curve  closes  there 


make_curve{) 

copy_curve() 

deletejcurvei) 

Curvei) 

print_curve() 

is_closed_curve 


Allocates  a  CURVE,  adds  to  INTERFACE 
Copies  a  CLUVE  to  current  INTERFACE 
Deletes  a  CURVZ  from  INTERFACE 
Create  and  add  CL'RVT  of  n  given  points 
Outputs  a  Formaned  CL'R\'E 
Checks  if  a  CLTIVT  is  closed 
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Pointi) 

copy_point() 

separationO 


Allocates  a  POINT  in  ourent  INTERFACE 
Copies  a  POINT  to  current  INTERFACE 
Computes  separation  of  two  POINTS 


BondO 

copyJbondO 

bondJengthO 


Allocates  a  BOND  in  current  INTERFACE 
Copies  a  BOND  to  current  INTERFACE 
Gives  the  length  of  a  BOND 


insert_point_in_bond() 
delete  _start_of_bond() 


Inserts  a  POINT  in  middle  of  a  BOND 
Deletes  POINT  at  start  of  BOND 


next_node() 
next_curve() 
nextjbondO 
next_point() 


Finds  Eest  >:ODE  on  an  INTERFACE 
Finds  Ltxt  CURVE  on  an  INTERFACE 
Finds  ne3rt'BOhT)  on  an  INTERFACE 
Finds  next  POINT  on  an  INTERFACE 


split j:urvei) 
join_cur\'es( ) 
intersectionsO 


Splits  a  CURVE  at  a  POINT  into  2  CL^RNTS 
Joins  CURVES  at  a  common  NODE 
Determines  Intersections  of  INTERFACE 


setJbaundaryO 


Finds,  Adds  Boundary  CLUVES  to  INTERFACE 


componentf)  Finds  COMPONENT  of  a  point  x,y 

max_componenti)  Finds  largest  COMPONENT  value 

min_componenti}  Finds  smallest  COMPONENT  value 

number_ofJabeled_componentsi)  Finds  Range  of  COMPONENT  values 


In  addition,  a  set  of  routines  is  available  that  allows  the  user  to  impose  a 
regular  rectangular  grid  on  some  rectangle  containing  or  contained  in  the 
region  where  the  interface  is  defined.  These  routines  allow  efficient  deter- 
mination of  the  topological  components  which  are  within  any  grid  block. 
They  assume  that  only  positive  COMPOr^TENT  values  are  used,  and  use  the 
negative  value  ONFRONT  to  indicate  grid  blocks  intersected  by  the  INTER- 
FACE.   A  related  set  of  routines  provide  information  about  the  BONDS  and 
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CURVES  in  a  grid  block. 


Operations  with  Superimposed  Grids: 


set_interface_grid() 

CompO 

number  j}f_comps() 

compsjnjblocki) 

is_component() 

tota  l_of_n  umber  _of_comp  s() 

number  jofJbondsO 

boruis_in_biocki] 

number _ofj:ur\es() 

curves_in_blockf) 


Establishes  a  regiJar  rectanguJar  grid 
Unique  grid  block  CONtPONENT  or  ONFRONT 
Number  of  CONfPONESTS  in  a  grid  block 
List  of  CONfPONENTS  in  a  grid  block 
Checks  if  a  COMPONENT  meets  grid  block 
Totals  over  ONTRONT  grid  blocks 
Number  of  BONDS  in  grid  block 
List  of  BONTO  in  ONTRONT  grid  block 
Number  of  CI,'R\'ES:in  grid  block 
List  of  Gi;R>I§dn-)Qf.TBONT  grid  block 


"tie.  fc  li.  S-'T.':  &  e:Mq' 

:  --L-30D  £  if   c3"'".."   CTj' 


-■  .-:r.'-:tj.!C>d  ^.n,-  .c. 
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Figure  1 :  The  fluid  flaw  in  a  two  dimensional  horizontal  oil  reservoir.  An  oil- 
water  interface  is  shown  at  a  selected  time  step,  resulting  from  the  injection  of 
water  into  five  wells  near  the  bottom  of  the  figure  and  the  production  of  oil  from 
five  wells  near  the  top  of  the  figure.  This  run  was  taken  from  unpublished  Joint 
work  of  the  authors  with  B.  Lindquist  and  G.  Tryggvason,  whom  we  thank  for 
permission  to  include  this  figure. 
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